# file = xcun.r  ssw 11/27/86  extend cunningham tables

define(STDOUT,6)
define(MPDIM,80)
define(MPCDIM,400)

implicit undefined (a-z)
common /xxx/col,b,n,pm,start,dcta,dta,won,pplim,tdlim,ppart
integer i,j,dcta,dctb,dta(50),dtb(50),iabs,jacobi,mod,col,idec
integer tdlim,ppart(1600,40),LMlim
logical islm

integer b,n,pm,del,LM,LMmod,LMres,start,limit,won(2),pplim
open(1,file='appc')
10 read(5,100,end=99) b,n,pm,del,LM,LMmod,LMres,start,limit,LMlim
100	format(10i7)
write(6,100) b,n,pm,del,LM,LMmod,LMres,start,limit,LMlim
# b=base, n=first entry (usually 1), pm=+-1, del=increment (1,2 or 4)
# LM=0 if table has no LM's, LM=1 if table has LM's
# n has LM iff  n mod LMmod = LMres
# do trial division up to tdlim iff start <= n <= LMlim
# do only LM's for limit < n <= LMlim

tdlim=100000
pplim=1+LMlim/3 # store pp iff n <= pplim to save space
won(1)=2 ; won(2)=1 # won = 1

for (; n <= LMlim; n=n+del) { # begin main loop

islm = (LM == 1) & (mod(n,LMmod) == LMres) # true if n has LM's

if (n > limit & ! islm) next # do only LM after limit
# form a list of divisors of n whose codivisors are odd
dcta=0 #count forward list
dctb=0 #count reverse list
# omit d=n/n if n is even
if (and(n,1) == 1 & n > 1)
	{
	dcta=1
	dta(1)=1
	}
do i=2,n
	{
	j=n/i # j = cofactor of i
	if (j < i) break
	if (n != i*j) next
	if (and(i,1) == 1)
		{
		dctb=dctb+1
		dtb(dctb)=j
		}
	if (j <= i) break
	if (and(j,1) == 1)
		{
		dcta=dcta+1
		dta(dcta)=i
		}
	}
# now copy the reverse list onto the forward list
if (dctb > 0)
	{
	do i=1,dctb
		{
		dcta=dcta+1
		dta(dcta)=dtb(dctb+1-i)
		}
	}

call putdec(STDOUT,n,5)		# write n
call putc(' ') ; call putc(' ') # and a couple of blanks
col=7				# column count

if (islm) # if n has LM's
	{
	dctb=0 # count entries in header line
	do i=1,dcta
		{
		j=jacobi(b,n/dta(i))
		if (j == 0) 
			{
			dctb=dctb+1
			if (dctb == 1) call putc('(')
			else call putc(',')
			call putdec(STDOUT,dta(i),1)
			if (mod(dta(i),LMmod) == LMres)
				{
				call putc('L'); call putc(',')
				call putdec(STDOUT,dta(i),1); call putc('M')
				}
			}
		dta(i) = j*dta(i) # either set sign or kill it
		}
	call putc(')'); call putc(' '); call putc('L')
	call putc('.'); call putc('M'); call putc(char(10)) # write ") L.M\n"

	do i=1,5
		call putc(' ')
	call putc('L'); call putc(' ')
	j=0 # count non-zero entries in dta
	do i=1,dcta
		{
		if (dta(i) != 0)
			{
			j=j+1
			if (j == 1) call putc('(')
			else call putc(',')
			col=col+idec(iabs(dta(i)))+2
			if (dta(i) > 0) call putc('L')
			else call putc('M')
			}
		}
	if (j > 0) {call putc(')'); col=col+1}
	col=col+1
	call putc(' ')
	call proc(1)

	do i=1,5
		call putc(' ')
	call putc('M'); call putc(' ')
	col=7
	j=0 # count non-zero entries in dta
	do i=1,dcta
		{
		if (dta(i) != 0)
			{
			j=j+1
			if (j == 1) call putc('(')
			else call putc(',')
			col=col+idec(iabs(dta(i)))+2
			if (dta(i) > 0) call putc('M')
			else call putc('L')
			}
		}
	if (j > 0) {call putc(')'); col=col+1}
	col=col+1
	call putc(' ')
	call proc(-1)
	}
else # else n does not have LM, so just dump the ()
	{
	if (dcta > 0)
		{
		do i=1,dcta
			{
			if (i == 1) call putc('(')
			else call putc(',')
			col=col+idec(dta(i))+1
			}
		call putc(')')
		col=col+1
		}
	call putc(' ')
	col=col+1
	call proc(0)
}

} # end main loop
goto 10
99 stop
end

integer function idec(n) # output n and count columns used
implicit undefined (a-z)
character ch(79)
integer itoc,i,n,nd
nd=itoc(n,ch,79)
do i=1,nd
     call putc(ch(i))
idec=nd
return
end

subroutine apcdec(k)
# output k to appc in i4 format
integer k,i,itoc,j
character lildig(10)

i=itoc(k,lildig,10)
do j=i,3
	call fputc(1,' ')
do j=1,i
	call fputc(1,lildig(j))
return
end

subroutine proc(m)
implicit undefined (a-z)
common /xxx/col,b,n,pm,start,dcta,dta,won,pplim,tdlim,ppart
integer m,col,b,n,pm,start,dcta,dta(50),won(2),pplim,ppart(1600,40)
integer t,i,j,k,ppt(MPDIM),tmp(MPDIM),quot(MPDIM),rem(MPDIM)
integer div1,gcd1,tdlim,idec,mpitoc
character bigdig(MPCDIM)

if (m == 0) # not LM here
	{
	call mpow(b,n,ppt)
	if (pm == 1) call add(ppt,won,ppt)
	else call sub(ppt,won,ppt) # ppt = b^n +- 1
	if (dcta > 0)
		{
		do i=1,dcta
			{
			k=ppart(dta(i),1)
			do j=1,k
				tmp(j)=ppart(dta(i),j)
			call div(ppt,tmp,quot,rem)
			if (rem(1)>2 | rem(2)>0) # if rem != 0
				{
				write(6,101) b,n,m,i,dta(i)
101				format('error dividing cofactor',5i12)
				stop
				}
			call mpcopy(quot,ppt)
			}
		}
	if (n <= pplim)
		{
		do j=1,ppt(1)
			ppart(n,j)=ppt(j)
		}
	}
else if (m == 1) # L here
	{
	call aurif(b,n,m,ppt)
	do i=1,dcta
		{
		t=iabs(dta(i))
		if (t == 0) next
		if (dta(i) < 0) t=t+800
		k=ppart(t,1)
		do j=1,k
			tmp(j)=ppart(t,j)
		call div(ppt,tmp,quot,rem)
		if (rem(1)>2 | rem(2)>0) # if rem != 0
			{
			write(6,101) b,n,m,i,dta(i)
			stop
			}
		call mpcopy(quot,ppt)
		}
	if (n <= pplim)
		{
		do j=1,ppt(1)
			ppart(n,j)=ppt(j)
		}
	}
else # M here
	{
	call aurif(b,n,m,ppt)
	do i=1,dcta
		{
		t=iabs(dta(i))
		if (t == 0) next
		if (dta(i) > 0) t=t+800
		k=ppart(t,1)
		do j=1,k
			tmp(j)=ppart(t,j)
		call div(ppt,tmp,quot,rem)
		if (rem(1)>2 | rem(2)>0) # if rem != 0
			{
			write(6,101) b,n,m,i,dta(i)
			stop
			}
		call mpcopy(quot,ppt)
		}
	if (n <= pplim)
		{
		do j=1,ppt(1)
			ppart(n+800,j)=ppt(j)
		}
	}

# check for intrinsic factor (at most one prime)
t=gcd1(div1(ppt,n,tmp),n) # t = gcd(ppt,n)
if (t > 1) # if we have an intrinsic factor
	{
	col=col+idec(t)+2
	call putc('*') ; call putc('.')
	j=div1(ppt,t,ppt) # divide t out of ppt; j is dummy
	}

if (n < start) # if n is too small to bother with, quick finish
	{
	if (n < 40) call mptdec(STDOUT,ppt,1)
	else {call putc('o'); call putc('m'); call putc('i'); call putc('t')}
	}
else	# else do some trial division and output to appc
	{
	if (n==1)
		while (and(ppt(2),1)==0)
			{
			col=col+idec(2)+1
			call putc('.')
			ppt(2)=ppt(2)/2
			}
	k=n+n
	for (t=k+1; t<tdlim; t=t+k)
		{
		while(div1(ppt,t,tmp)==0)
			{
			col=col+idec(t)+1
			call putc('.')
			call mpcopy(tmp,ppt)
			}
		if (ppt(1)==2 & ppt(2)==1) break
		}
	do i=col,73	# space over to column 74
		call putc(' ')
	call putc('C')
	k=mpitoc(ppt,bigdig,MPCDIM)
	i=idec(k)	# write size of ppt to table
	call apcdec(k)	# write size of ppt to appc
	call apcdec(b)	# write base to appc
	call fputc(1,',')
	call apcdec(n)	# write n to appc
	if (m==0)	# write +-LM to appc
		if (pm==1) call fputc(1,'+')
		else call fputc(1,'-')
	else
		if (m==1) call fputc(1,'L')
		else call fputc(1,'M')
	call fputc(1,' ')
	call fputc(1,' ')
	do j=1,k
		{
		if (j > 1 & mod(j,100) == 1)
			{
			call fputc(1,char(10))
			call apcdec(k)
			do i=1,12
				call fputc(1,' ')
			}
		call fputc(1,bigdig(j))
		}
	call fputc(1,char(10))
	}

#if (ppt(1)>2 | ppt(2)!=1) call mptdec(STDOUT,ppt,1)
#else call putc('X')
call putc(char(10))
return
end

subroutine mpow(b,n,pow)
implicit undefined (a-z)
integer b,n,pow(MPDIM) # compute pow := b^n
integer e,i,j,bits(30),tmp(MPDIM)

pow(1)=2
pow(2)=b
e=n
i=0
do j=1,30
	{
	i=i+1
	bits(i)=and(e,1) # get low order bit
	e=e/2
	if (e==0) break
	}

for (i=i-1; i>0; i=i-1)
	{
	call mul(pow,pow,tmp)
	if (bits(i) != 0) call mul1(tmp,b,tmp)
	call mpcopy(tmp,pow)
	}
return
end

subroutine aurif(b,n,m,ppt)
integer b,n,m,ppt(MPDIM),w(MPDIM),x(MPDIM),y(MPDIM),z(MPDIM),h,k

h=n/b
k=(h+1)/2

switch (b) {

case 3:
# aurif for b = 3 ; n = 3h and h = 2k-1
# L = 3^h - 3^k + 1  for m=1
# M = 3^h + 3^k + 1  for m=-1

call mpow(3,h,x) # x = 3^h
call mpow(3,k,y) # y = 3^k
x(2)=x(2)+1      # x = x+1
if (m==1) call sub(x,y,ppt)
else call add(x,y,ppt)

case 5,6:
# aurif for b = 5 or 6 ; n = b*h and h = 2k-1
# L = x^2 - x*b^k + b^h for m=1
# M = x^2 + x*b^k + b^h for m=-1
# x = b^h + 1

call mpow(b,h,z)  # z = b^h
call mpcopy(z,x)  # x = z
x(2)=x(2)+1       # x = b^h + 1
call mul(x,x,y)   # y = x^2
call add(y,z,y)   # y = x^2 + b^h
call mpow(b,k,z)  # z = b^k
call mul(x,z,w)   # w = x*b^k
if (m==1) call sub(y,w,ppt)
else call add(y,w,ppt)

case 7:
# aurif for b=7 ; n = 7h and h = 2k-1
# L = w^3 - z for m=1
# M = w^3 + z for m=-1
# w = 7^h + 1
# z = (w^2 - 7^h)*7^k 

call mpow(7,h,w)  # w = 7^h
w(2)=w(2)+1       # w = 7^h + 1
call mul(w,w,x)   # x = w^2
call mul(w,x,ppt) # ppt = w^3
w(2)=w(2)-1       # w = 7^h
call sub(x,w,x)   # x = w^2 - 7^h
call mpow(7,k,y)  # y = 7^k
call mul(y,x,z)
if (m==1) call sub(ppt,z,ppt)
else call add(ppt,z,ppt)

case 10:
# aurif for b=10 ; n = 10h and h = 2k-1
# L = y - x (for m=1) and M = y + x (for m=-1)
# y = 10^4h + 5*10^3h + 7*10^2h + 5*10^h + 1
# x = (10^3h + 2*(10^2h + 10^h) + 1)*10^k

call mpow(10,h,z)  # z = 10^h
call mul1(z,5,x)   # x = 5*10^h
call mul(z,z,y)    # y = 10^2h
call mul(y,y,w)    # w = 10^4h
call add(y,z,z)    # z = 10^2h + 10^h
call add(z,z,z)    # z = 2*(10^2h + 10^h)
call mul1(y,7,y)   # y = 7*10^2h
call add(x,y,y)    # y = 7*10^2h + 5*10^h
call add(w,y,y)    # y = 10^4h + 7*10^2h + 5*10^h
z(2)=z(2)+1	   # z = z + 1
y(2)=y(2)+1	   # y = y + 1
call mpow(10,3*h,w) # w = 10^3h
call add(w,z,z)    # z = 10^3h + 2*(10^2h + 10^h) + 1
call mul1(w,5,w)   # w = 5*10^3h
call add(w,y,y)    # y = final
call mpow(10,k,w)  # w = 10^k
call mul(w,z,x)    # x = final
if (m==1) call sub(y,x,ppt)
else call add(y,x,ppt)

case 11:
# aurif for b=11 ; n = 11h and h = 2k-1
# L = z - y (for m=1) and M = z + y (for m=-1)
# z = 11^5h + 5*11^4h - 11^3h - 11^2h + 5*11^h + 1
# y = (11^4h + 11^3h - 11^2h + 11^h + 1)*11^k

call mpow(11,h,z)  # z = 11^h
call mul(z,z,y)    # y = 11^2h
call mul(y,y,x)    # x = 11^4h
call add(x,z,x)    # x = 11^4h + 11^h
call mul1(x,5,z)   # z = 5*11^4h + 5*11^h
call sub(x,y,x)    # x = 11^4h - 11^2h + 11^h
call sub(z,y,z)    # z = 5*11^4h - 11^2h + 5*11^h
x(2)=x(2)+1	   # x = x + 1
z(2)=z(2)+1	   # z = z + 1
call mpow(11,3*h,w) # w = 11^3h
call sub(z,w,z)    # z = z - 11^3h
call add(x,w,x)    # x = x + 11^3h
call mpow(11,5*h,w) # w = 11^5h
call add(z,w,z)    # z = final
call mpow(11,k,w)  # w = 11^k
call mul(w,x,y)    # y = final
if (m==1) call sub(z,y,ppt)
else call add(z,y,ppt)

case 12:
# aurif for b=12 ; NOTE n = 3h and h = 2k-1 NOTE
h=n/3
k=(h+1)/2
# L = 12^h - 2^h*3^k + 1 for m=1
# M = 12^h + 2^h*3*k + 1 for m=-1

call mpow(2,h,x)  # x = 2^h
call mpow(3,k,y)  # y = 3^k
call mul(x,y,z)   # z = 2^h*3^k
call mpow(12,h,y)  # y = 12^h
y(2)=y(2)+1	  # y = y + 1
if (m==1) call sub(y,z,ppt)
else call add(y,z,ppt)

default:
write(6,100) b
100	format('aurif not defined for b =',i10)
stop

}
return
end

integer function gcd1(u,v)
implicit undefined (a-z)
# given u,v, find gcd(u,v)
integer u,v,u3,v3,t,mod

u3=u
v3=v

while (v3!=0)
	{
	t=mod(u3,v3)
	u3=v3
	v3=t
	}
gcd1=u3
return
end
